Numerical Studies of Density-driven Flow in Co2 Storage in Saline Aquifers

Simulations are routinely used to study the process of carbon dioxide (CO2) sequestration in saline aquifers. In this paper, we look at some numerical aspects of the accurate modeling and simulation of the dissolution-diffusion-convection process. We perform convergence studies with respect to solver tolerances, grid resolutions, fluctuation strength, and domain size. We show that stringent tolerances and grid resolutions are needed to accurately predict onset time. Domain size must be sufficiently large to contain at least 2 extended fingers to accurately predict the long-term stabilized mass flux of CO2; otherwise, finite domain effects will adversely change the flow behavior of the system we are modeling. INTRODUCTION Carbon dioxide (CO2) sequestration involves injecting CO2 into a saline aquifer. While the primary mechanism of securing the CO2 relies on a leak-proof formation, secondary geochemical mechanisms may play a significant role, especially in a geological time frame. At long time, an immiscible CO2 gas layer will form on top of the brine in the rock formation. Under ambient temperature and pressure conditions in a typical aquifer, CO2 will dissolve into the brine and increase the density of the brine at the interface of the layers by 0.1–1%, depending on the salinity of the brine (Pruess and Zhang, 2008). Due to gravitational instability and the heterogeneity in the rock properties of the aquifer, CO2-rich brine fingers will form, leading to convective flow that transports these CO2-rich brines downward, while driving brine with low CO2 concentration upwards. This then accelerates the rate at which CO2 is dissolved and provides a more secure mechanism by which CO2 can be stored. This dissolution-diffusion-convection process has been analyzed in a number of studies. In Ennis-King and Paterson (2003) and Riaz et al. (2006), linear stability analyses yield useful relations for the onset time for convection, dominant wavelength for growth of convective fingers, and the growth rates of these fingers. Numerical simulations were also performed to further elucidate the process and to validate the linear stability analyses. For example, Riaz et al. (2006) performed numerical simulation of a singlephase two-component model with the Boussinesq assumption and demonstrated that the simulation results are consistent with their analysis. Pruess and Zhang (2008) examined long-term behavior of the CO2 flux, in addition to the onset of convection. Their simulation uses a full compressible model with very accurate equations of state. In this work, we examine some of the numerical aspects of studying the dissolution-diffusionconvective process of CO2 through simulation. We use a second-order accurate adaptive method that is described in the next section. Specifically, we examine how simulation parameters, such as solver tolerances, grid resolution, strength of perturbations, and domain size, affect our solution. The results are compared to those obtained through TOUGH2-MP, a parallelized version of the general-purpose simulator TOUGH2/ECO2N (Zhang et al., 2008; Pruess, 2004; Pruess and Spycher, 2007).